{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cd172d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Created on Mon May 16 19:00:32 2022\n",
    "@author: Ajit Johnson Nirmal\n",
    "SCIMAP tutorial May 2022\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7a1c1cc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load packages\n",
    "import scimap as sm\n",
    "import scanpy as sc\n",
    "import pandas as pd\n",
    "import anndata as ad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e2f30407",
   "metadata": {},
   "outputs": [],
   "source": [
    "common_path = \"/Users/aj/Dropbox (Partners HealthCare)/conferences/scimap_tutorial/may_2022_tutorial/\"\n",
    "#common_path = \"C:/Users/ajn16/Dropbox (Partners HealthCare)/conferences/scimap_tutorial/may_2022_tutorial/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecd0fb93",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data\n",
    "#adata = sm.pp.mcmicro_to_scimap (image_path= str(common_path) + 'exemplar_001/quantification/unmicst-exemplar-001_cell.csv')\n",
    "#manual_gate = pd.read_csv(str(common_path) + 'manual_gates.csv')\n",
    "#adata = sm.pp.rescale (adata, gate=manual_gate)\n",
    "#phenotype = pd.read_csv(str(common_path) + 'phenotype_workflow.csv')\n",
    "#adata = sm.tl.phenotype_cells (adata, phenotype=phenotype, label=\"phenotype\") \n",
    "# add user defined ROI's before proceeding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "360428e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load saved anndata object\n",
    "adata = ad.read(str(common_path) + 'may2022_tutorial.h5ad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "59758fee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AnnData object with n_obs × n_vars = 11170 × 9\n",
       "    obs: 'X_centroid', 'Y_centroid', 'Area', 'MajorAxisLength', 'MinorAxisLength', 'Eccentricity', 'Solidity', 'Extent', 'Orientation', 'imageid', 'phenotype', 'index_info', 'ROI', 'ROI_individual'\n",
       "    uns: 'all_markers', 'dendrogram_phenotype'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "adata"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad48620c",
   "metadata": {},
   "source": [
    "### Investigate cell-type composition within the ROI's"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5c8ca0f8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAE9CAYAAACFlCHjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABHq0lEQVR4nO3deXiU1fn/8fckk4SwhVV2wY1bQIWKiKIoboCi1AVFxCJYwVbLV/tTbFFxq1apolate10oKi61irKpqAWlKLjgAt4oBhTZgmxhyTKZ/P6YCYQYhgCBeUI+r+viInOeZT4Twsydc87znFBxcTEiIiIi25OS7AAiIiISbCoWREREJCEVCyIiIpKQigURERFJKJzsADtiZhlAF2AZUJTkOCIiVUUq0AyY7e75yQ4jVVvgiwVihcKMZIcQEamiugMfJDuEVG1VoVhYBvDcc8/RtGnTZGcREakSli9fzsCBAyH+HlraJ598kh4Oh58AjifWAyFSBHwQiUSGdu7cuaDsxqpQLBQBNG3alJYtWyY7i4hIVfOL4duUlJTf161b97jWrVuvTUlJ0c12hGg0Glq8ePHxa9eu/T3w97LbNcFRRKSaSU1NHdK8efONKhSkREpKSnHz5s03pKamDi53+17OIyIiSVZcXJyVnp5emOwcEizp6emFxcXFWeVtU7EgIlL9hEKhULIzSMDEfybKrQuqwpwFERHZgyKRSMdwOFzpnweRSCQSDofnVvZ5Ze9TsSAiUs2Fw+Hw448/XunnHTZsWIU+Y7Kzs9MHDRpkM2bM+LJ0u5l1dvdPyjvm/fffr/PQQw81f+WVV7wyskpiFfqHNLOLgZHxh5Pd/Voz6wQ8AWQB04HfuXvEzPYHxgH7AQ4MdPcNZlYPeA44EMgBLnD35ZX5YkRERKTy7bBYMLOawANAW2At8KGZnQrcD1zm7rPM7J/AUOAR4GHgYXcfb2ajgFHAn4DbgRnu3sfMfkPs0oz+lf+SRERkXzF27NiGH3zwQVZubm7qsmXLMrp06bL+7rvv/qH0Pv/4xz/2e++99+o/++yz315yySWHtG/ffuPcuXNrr1u3LvynP/3ph9NPP3398uXLwyNGjGizcuXK9NTU1OL/+7//++nII4/c1Ldv3/azZs36AuDYY4894o9//OOPF1xwwZoxY8Y0TUlJKc7Ly0tduXJl2o8//lhjxYoV6X379s0ZMWJEtftFtyITHFPj+9UC0uJ/CoFMd58V3+cZ4HwzSwNOAF4p3R7/ug+xngWAF4DT4/uLiIhs19dff13rscceWzhx4sR5M2fOrPfFF19klmx79tlnG7777rv1n3nmmW9r1aoVBYhEIqHXX3/9m2uuuebHBx98sAXAqFGj9j/66KNzp06dOu/BBx9ceMstt7SJRqPst99+BV9++WWN+fPn1ygqKgrNnj27DsDMmTOzevbsuQ7gu+++y3zuuecWvPLKK/PHjRvXbM2aNdXuRlY7LBbcPZdY78A3wE/AIqCAbe8KtgxoCTQC1rt7pEw7QPOSY+Lb1wONd/sViEiVF4lEdrzTXhS0PPu6lJRffhRFo9GS2fl06NBhY926daO1atWKNm3aNH/16tWpAIsWLaoxevTo1gMHDlxRu3btaMmxJ5xwwrr4cZtzc3PDAJ999lmd3/zmN6sADjrooIJ27dptnD17dq3jjjtu3fTp0+vOmDGjzgUXXLBi7ty5tdeuXZu6Zs2atA4dOuQBdO7cOTcjI6O4SZMmkTp16kTWrVtX7YqFigxDHAFcCrQG1hGbj9CznF2jQHnX4pT8AybaJiLVWDgcZk9MsNtVw4YNS3aEaqV+/fqRjRs3bvMBvHLlynDt2rWLADIyMrZ8VsQLiBBAZmZm9Oabb150991379+zZ8/1JQVDjRo1ikvtC0BxcfE2n0HFxcVEIpHQKaecsu6BBx5onp6eHr3mmmuWvvPOOw1eeumlBkcfffS6kn0zMjK23LwqFApRXFz97mVVkWGIXsA0d18ZX7nsGaAHUHqhhmbAUmITF+uaWWqZdoj1SjQFMLMwUBf4eTfzi4hIFVe3bt1oixYt8l999dV6JW3PPvts4yOPPHJ9ouOaNGlScNZZZ63r1KlT7l133dU80b6dOnVa/69//asRwMKFC9O//vrr2l27dt145JFHblqyZEmNH3/8sUa7du3yjjzyyNynn3662cknn7wu0fmqm4pcDTEX+JuZ1QI2AWcB/wX6mdlx7v4hMIjYVRKFZjaD2MTF50va4+eZFH/81/j2Ge6uO4iJiCRZJBKJVPQyx509b0Vv33DPPfd8P2rUqNaPPfZY80gkEjrooIM233HHHT9Mnjy53DsKlnbTTTf92KdPn8POOeec1dvb57bbbvvxz3/+c+uePXs2DIVCjBo1anHz5s0LAQ4//PDczZs3pwIce+yx6ydOnNjoxBNPzK3gy6wWQhXpTjGzPwFDiE1s/Bi4EjBil07WAT4Dhrh7vpm1Bp4ldunkD8AAd19jZg2I9UocROyqioHuvqgCz90GyJ42bZoWkhLZh2kYonItWbKEU045BeCAsu+1c+fOXdSxY8dVSQkmgTZ37txGHTt2bFO2vUIln7uPBkaXPSdwdDn7LiY2TFG2fTXQtyLPJyIiIsGhtSFEREQkIRULIiIikpCKBREREUlIxYKIiIgkpGJBREREEtIS1SIi1VxxYVHHUFpqpX8eFBcWRUJpqXMT7TNixIj9v/zyy9qRSCS0bNmyjFatWuUBDBgwYMUll1yiG/cFhIoFEZFqLpSWGs677f1KP2+Nm3rs8DOmZAXJ7Ozs9EGDBtmUKVPmVXoQ2W0ahhARkcDp3r374dnZ2ekA77//fp1+/foZQL9+/ez6669v1bNnz/Y9evQ4bPLkyXUHDhx4yHHHHXf4gw8+uB/Axo0bU37/+98fcNppp3Xo1atX+3HjxjWE2HLX559/ftuePXu2v+2221ok79VVPepZEBGRKuett96ad9dddzUbPXr0/hMnTpyXk5MTPvfcc9sPHz585ejRo5tnZWVF3n777a9zcnLC5513XrvDDz98E0BOTk7622+//VVaWlqyX0KVop4FERGpUnr06LEOoHnz5gXt27ffWKtWrWibNm0KSlau/OSTT+pceOGFqwAaN24cOf7449d++OGHdQDatm27SYXCzlOxICIigVSydlFhYeE2y0unp6dvWdQoNTX1FwsclbfmUVFRUQi2Xe5aKk7FgoiIBE5WVlZk3rx5mQBvv/12vZ05tnPnzrnjx49vBJCTkxOeMWNGvW7dumkVyd2gOQsiItVccWFRpCJXLuzKeUNpqbt07BVXXLF09OjRrR599NHmXbt2Xbczx1577bVLR44c2bpnz57to9FoaMiQIcs6d+686euvv87cpTBSsSWqk0lLVItUD1qiunJpiWrZFdtbolrDECIiIpKQigURERFJSMWCiIiIJLTDCS1mdhnwh1JNBwD/Al4D7gUygRfd/cb4/p2AJ4AsYDrwO3ePmNn+wDhgP8CBge6+odJeiYiIiOwRO+xZcPcn3b2Tu3cCBgIrgdHAU8CvgXZAFzM7PX7IOGC4u7cFQsDQePvDwMPufigwBxhVmS9ERERE9oydHYZ4BLgeOBD41t2z3T1CrEA438xaA5nuPiu+/zPx9jTgBOCV0u27mV1ERET2ggpfV2tmpxIrBF42swHAslKblwEtgebbaW8ErI8XFqXbRUQkyYoK8jumpmdU+n0WigryI6npGQmXqJaqYWd+OC4nNkcBYsMLZUV3oV1ERJIsNT0j/EL/dpV+3gEvzt/hZ8yIESP2//LLL2tHIpHQsmXLMlq1apUHMGDAgBWXXHLJz5UeSnZJhYoFM0sHTgQGx5t+ApqW2qUZsDRBew5Q18xS3b2oVLuIiFRjd9999w8A2dnZ6YMGDbIpU6bMS3Ym+aWK9iwcASxw943xxx8BZmYHA9nARcBT7r7YzPLM7Dh3/xAYBEx290IzmwH0B54vaa/UVyIiIvuM7t27H96uXbuN3333Xc0XXnjhm7fffjvr+eef3y8ajYbMbONdd931Q2ZmZvHLL79c/5FHHmleo0aNaNu2bTcVFRWFHnzwwUXJzr+vqegExwOBJSUP3D2PWC/Dv4F5wDdsnbw4ELjPzOYDtYAH4u1XAMPMbB7QHbhxd8OLiMi+6/jjj1/37rvvfpWTkxN+9dVXG/373//+ZsqUKfMaNGgQefDBB5uuXLkyPGbMmFZjx45d8MYbb8xfv3691jvaQyr0jXX3l4CXyrRNAzqWs+9c4Ohy2hcDPXYppYiIVDtHHXXURoAPPvigzk8//VTjnHPOaQcQiURCbdu23fThhx/W7tChw8aWLVsWApx99tk/7+wKlVIxqsJEJOkikUigFm+KRCKEw3p7TLbMzMwoQFFRUeikk05afdddd/0IkJubmxKJREIzZsyoE41qrvzeoP8NIpJ04XCYvNveT3aMLWrc1CPZEaSU4447LvfKK69sumLFimWNGzeO/OlPf2rdqlWr/N/+9rcr77rrrv2XLl2a1rRp08KJEyfWT09PD/ZSylWUigURkWquqCA/UpHLHHflvKnpGbt9nk6dOm2+7LLLll588cVWXFzMIYccsunqq69elpmZWTxixIgfLrnkkrbp6enRpk2b5tetW7eoEqJLGSoWRESquT1146SdKRQOOOCAghkzZnxZ8rj01wBDhgxZNWTIkFWl21atWpU6f/78mlOmTPk6NTWVP//5z61at26dv9vB5RdULIiISJXUoEGDotzc3NRevXp1SE1NLW7btu2mwYMH5yQ7175IxYKIiFRJKSkp3HnnnT8mO0d1sLMLSYmIiEg1o2JBREREElKxICIiIgmpWBAREZGENMFRRKTay+8IGXvg8yA/Aju+LDM3Nzfltttuazl79uy6GRkZ0Vq1ahVdeeWVS0855ZRcgH/+85+NateuHe3fv//q4cOHt+nSpUvuoEGDdnv56rFjxzacPXt2nX1t4amSFTxnzJjxZWV9v1QsiIhUexlhOGoPnHfODj9jotEol1566cGHHHLI5qlTp36dkZFR/Omnn2ZeccUVh6Smpmb36NEj9/PPP6/dpUuX3D0QUCpIxYKIiCTN9OnT66xYsSLjxRdfXJCSEhsZP/LIIzdfeumlyx5++OFmkUiEmTNn1vv000/rNGnSpDB+TNbLL7+835o1a8JDhgxZ9tvf/nZVbm5uysiRI/f//vvvM6PRaGjIkCHL+/fvv3rs2LEN33jjjYbr1q0LH3/88etuuummn8rLMXz48DY1atSIfvHFF7U3btyYeu211/44YcKEhgsXLszs3r372ttvv33J2LFjG06fPj1r1apV6Tk5OWkXXnjhymXLlqV/8sknderWrRsZO3bst8uXL08r+a0e4M4772wOMHLkyKXHHHPMESeeeOKaL774ok5qamrxAw88sPDAAw8s+Pjjj2veeeedrfLz81OysrIid9xxx+IDDzywoHS+zz77LPOmm25qnZ+fn1KnTp2i+++///tWrVoVjhkzpum0adPqFxUVhbp27br+lltuWVLe61u3bl3KlVdeeeDq1avTAH7/+98vPeuss9ZV9N9JcxZERCRpPv/885pmtrGkUChx7LHH5i5YsKDWqaeemtutW7e1l19++dJevXqtBygoKEh5/fXX5z/yyCPfPvrooy0AxowZ06x9+/abJk2aNP/FF1/85sknn2y2cOHCdICcnJz0iRMnztteoVBi1apVaVOnTp13+eWXL7311lvb3HnnnYvfeOONeW+++WbjtWvXpgJ88803tZ599tkFzz33nD/00EMtTzzxxHVTp06dB/D2229nJTr/mjVr0rp165Y7efLkeZ06dcp96qmn9svPzw+NGjWqzf333589adKk+YMHD14xcuTINmWPve666w68/PLLl7311lvzevXqtfqJJ55oMmXKlLrz58+vOWHChPmTJk2al5OTkzZ+/PgG5T33hAkT6jdr1qxg0qRJ8++5557s2bNn10mUtSz1LIiISNKEQiGKiopCZdsLCgq2+8vsSSedtDYlJYUOHTrkrV+/Pgwwe/bsuvn5+SkTJkxoBJCXl5cyb968TIC2bdtuSktL22GW448/fh1Aq1atCtq0abO5SZMmEYA6depE1qxZkwpw2GGHbcjKyopmZWUVAJx88snrAZo2bVpQUlAk0rNnz3UAhxxyyOY5c+bUWbBgQcby5cszLr/88oNL9tm0adM258nJyQmvXr067cwzz1wHMGzYsByAUaNGtZw/f36tM888s33J96xp06YFxx577Iayz9u1a9cNDz/8cIvBgwennXDCCeuuueaapTv8hpSiYkFERJLmV7/61caXXnppv4KCglDpFSPnzJlTq23bthvLOyY1NbUYYndwLBGNRkOjR4/O7ty58yaA5cuXhxs0aFA0fvz4BhkZGRVaxzotLW3L86emlv+5X3qf+ONttodCIYqLt+4SiURC4XB4S0NmZmZxyX4QW367adOm+VOmTJkX358VK1Zsc9Kyz7l58+bQ0qVL04qKikIXXnjhyuHDh68AWLNmTWo4HC5etWrVLz7b27Ztmz916tSv3nrrraz33nsva9y4cU3eeeedr8v26GyPhiFERCRpTjjhhA1t2rTJu+GGG1oVFBSEAObMmVPz6aefbnbFFVcsg1hxUF7vQ2mdO3deP27cuMYAS5cuTevbt2+HH374IX3Pv4Jt1atXr2jDhg3hlStXhvPy8kL/+9//6iba/9BDD83Lzc0NT58+vTbA2LFjG1111VUHlj1no0aNCt5+++26AC+88ELDMWPGtDj22GPXT5o0qWFubm5KYWEhQ4cOPfg///lP/fKe59FHH208evTo5v369Vtz1113/bBu3bq0devW7bAnpESFehbM7CzgFqAWMNXdrzKzU4F7gUzgRXe/Mb5vJ+AJIAuYDvzO3SNmtj8wDtgPcGCgu/+iq0RERPa2/EhFrlzYtfPueOXJJ5544rs77rijRe/evTukpKQU16lTJ3LHHXdk9+jRIxfg2GOPXf+Pf/yjZaLlp0eMGLF05MiRrU877bQO0WiU4cOHLzn44IPzZ86cWbsSX9AO1atXr2jAgAHLzzvvvHaNGzcuaN++fbm9IyVq1KhRfM899yy8884797/jjjtCNWvWLLr77rsXld3v7rvvzr755pv3HzNmTMusrKzIfffdl928efPC+fPn1zz33HPbFRUVccwxx6y/+OKLf168ePEviqQBAwb8fOWVVx7Ys2fP9uFwuPi3v/3t0vr161d4Oe9Q6e6S8pjZgcAMoCuwAngX+CvwGHAi8CMwEbjf3Seb2VfAZe4+y8z+Ccxx90fM7E1gnLuPN7NRQG13/9OOAppZGyB72rRptGzZsqKvS0SqmLzb3k92hC1q3NQj2RF225IlSzjllFMADnD3RaW3zZ07d1HHjh1XlXugVGtz585t1LFjxzZl2ysyDHEOsZ6DJe5eCPQHNgHfunu2u0eI9Ricb2atgUx3nxU/9pl4expwAvBK6fbdeD0iIiKyl1Sk2+lgoMDMpgJNgTeAr4FlpfZZBrQEmm+nvRGwPl5YlG4XERGRgKtIsRAm1ivQA9gAvE6sZ6GsKFDeBJRE7SIiIhJwFRmGWA684+457r4ZeA04jVgvQ4lmwFLgp+205wB1zSy1TLuIiIgEXEWKhTeBXmZWL/5hfzqxuQdmZgfH2y4CJrv7YiDPzI6LHzso3l5IbJJk/9LtlflCREREZM/YYbHg7h8BfwM+AOYBi4FHgMHAv+Nt37B18uJA4D4zm0/sUssH4u1XAMPMbB7QHbix0l6FiIiI7DEVuq7W3Z8CnirTPA3oWM6+c4Gjy2lfTGzeg4iIBEikINIxnB6u9PssRAoikXB6eIdLVEvw6XbPIiLVXDg9HH7s1+Mq/byXv35xhT5jsrOz03v37n34Qw899O1pp522vqS9e/fuh48dO9YPOOCAgkTHlzV8+PA211xzzdI2bdrs1HHl2V6Gq6++uvXAgQNzNm7cmPrQQw81f+WVV3x3nyvIdLtnERFJunA4XHzbbbe1Xr9+/W5/Ln3++ed1dnTDwd11//33L+7SpUt5Vwbuk9SzICIiSdegQYPCo446av3NN9/c6r777ltcetuYMWOaTps2rX5RUVGoa9eu62+55ZYlixcvTh80aJDNmDHjS4A777yzOUB6enp09erVaUOHDj1k/Pjx35xzzjnt27Vrt/G7776r+cILL3wzZcqUrHHjxjUNhULFZrbpr3/96w916tSJdunSpeMxxxyzbsGCBTUzMzOj99133/clvQn33Xdf82+//bZmfn5+yl//+tfsY445ZmO/fv3sD3/4wzZX9T344INN3nzzzYYpKSm0a9du47333rvN66jK1LMgIiKBcMstt/w4e/bsuiULJgG88847WfPnz685YcKE+ZMmTZqXk5OTNn78+AbbO8c111yzvEGDBoVPPPHEt40aNSqC2NLT77777lcrVqxIe/rpp5u98MIL37z11lvzatSoEb377rubA6xfvz589NFH506dOnVer169Vt9yyy37l5zzoIMO2jx58uR5F1xwwYonn3yySXnPW1hYyNixY5u+8cYb89988815KSkpxUuWLNnxuthVhIoFEREJhKysrOhNN920uPRwxKxZs+rMnz+/1plnntm+T58+7d291rfffpu5M+c96qijNgLMnDmzdrdu3daVFBEDBw7MmTNnTh2A9PT04oEDB/4McNFFF636/PPP65Qcf8YZZ6wFMLO8devWldsjn5aWRvv27Tf07du33ejRo5tfcsklK1u2bFm4C9+GQNIwhIiIBEbPnj3XT548ef3NN9/cCiAajYYuvPDClcOHD18BsGbNmtRwOFz8888/h0vPS4hEIqFwOFzuRIXMzMxoyblKtxcXF1Oy9HUoFCpOSYn9/hyNRkMpKSlbzlVy3lAoRHFx8XaXyn766acXzpo1q9Z7772Xdfnll7f961//+n2PHj32idWV1bMgIiKBUjIcsXr16rSjjz46d9KkSQ1zc3NTCgsLGTp06MH/+c9/6terV69ow4YN4ZUrV4bz8vJC//vf/7YMXaSmphZHIpFffKgfd9xxuR9++GHWzz//nArw/PPPNz7yyCNzAfLz81PeeOONLIBx48Y1PProo9ftTOaVK1eGTz311A6HH3745htuuGHpUUcdtX7+/Pk1d+87ERzqWRARqeYiBZFIRS9z3NnzhtN3/rQlwxHDhw8/pHfv3mtzc3NTzz333HZFRUUcc8wx6y+++OKfU1JSGDBgwPLzzjuvXePGjQvat2+/seT4bt26rRs2bNgh//znPxeUPm/Hjh03Dx48ePmAAQMsEomEzGzT6NGjt0xCnDJlSv2HHnqoRcOGDQvvvffe7J3JvN9++0XOOeecnLPPPrtdRkZGtEmTJgUDBw7cZ5YBD+3py0t2l5m1AbKnTZtGy5ZaqFJkX5V32/vJjrBFjZt6JDvCbluyZAmnnHIKwAHuvqj0trlz5y7q2LHjPvNBVhnMrLO7f5LsHMk2d+7cRh07dmxTtl3DECIiIpKQigUREan21KuQmIoFERERSUjFgoiIiCSkYkFEREQSUrEgIiIiCalYEBGp5gqKCjoCnSv7T/y8CWVnZ6ebWefS60FAbGno7Ozs9Pfff79Ov379rKR9/fr1Kb/+9a8PvfHGG3frWvqxY8c2HD58eJvSz7U759vX6aZMIiLVXHpqerjva30q/bwTzp5Yoc+YkuWpu3bt+nXdunWj29svNzc3ZdCgQW1/9atf5d5yyy0/VV5S2ZEK/UOa2btAE6BkUYzLgYOAG4F04D53/0d831OBe4FM4EV3vzHe3gl4AsgCpgO/c/dIpb0SERGpkhItT11iw4YNKYMHDz7kqKOOWn/jjTcuLW+f8ePHN3jyySebhUIhDj300I1jxoxZnJ+fHxo5cuT+33//fWY0Gg0NGTJkef/+/VeXd/zcuXMzR40a1bqoqCiUnp4eHT169KK2bdvmV+Zrrap2OAxhZiHgUKCju3dy907AEuAO4HigIzDMzNqbWSbwFPBroB3QxcxOj59qHDDc3dsCIWBoZb8YERGpmspbnrpEfn5+ypAhQw5etGhR5u9///sV5R2/ZMmStHvvvbfV008/veDtt9/+uqioKDR58uSsMWPGNGvfvv2mSZMmzX/xxRe/efLJJ5stXLiw3CGHJ598ssmgQYNWTJw4cf6AAQNWzp49u1Zlv86qqiI9CwYUA5PNbD9ivQO5wLvuvhrAzF4B+gH/Bb519+x4+zjgfDObB2S6+6z4OZ8BbgUeqcTXIiJVVHFhUaBusVxcWEQoLTXZMaqV0stTd+3a9evS2xYsWFDzt7/97dLWrVvnjRgxos1TTz21sOzxH330Ue3DDjtsQ6tWrQoBHn744WyAxx9/vFl+fn7KhAkTGgHk5eWlzJs3r9wlrnv06LFu9OjR+8+YMaPuSSedtO7cc89dU/mvtGqqSLFQH5gG/J7Y0ML7wIvAslL7LAOOBpqX094yQbuICKG0VF7o3y7ZMbYY8OL8ZEeolsouT13i0EMP3Xjdddct27hxY8qZZ57Z/oknnmg8dOjQnNL7pKWlbbPQ0cqVK8MQW2569OjR2Z07d94EsHz58nCDBg2Kxo8f36Ds85933nlrunbtumHKlClZ//rXv5r897//zdresEh1s8NhCHf/n7sPcveN7r4K+CdwWzm7RokNL+xMu4iIyBall6cuaSspBGrVqhW98847s//xj3+0/Prrr2uUPq5z584b58+fX2vZsmXh+HlaTZw4sV7nzp3Xjxs3rjHA0qVL0/r27dvhhx9+KHcYYtiwYQfOmTOn1mWXXbbq//7v/35y931miendtcOeBTM7Hshw92nxphCwCGhaardmwFLgp51sFxGRJCsoKohU9MqFnT1veurOXZFYennq8rYfc8wxGy+44IIV/+///b8DX3vttfmZmZnFAC1atCi89tprfxw8eHDbaDQaOuywwzb85je/WbVx48aUkSNHtj7ttNM6RKNRhg8fvuTggw/OnzlzZu2y577iiiuW3XjjjW0ef/zxZqmpqVx33XU/7tIL3wftcIlqMzuTWE9CNyAN+AC4gtiExaOBjcBMYBjwBfAtcBKQDbwJPOXuL5vZV8Dl7v6hmT0BLHD3u3cUUEtUi1QPGoaoXFqiWnbFLi9R7e5vAhOBz4BPiH34fwjcALwHfA487+4fu3seMBj4NzAP+AZ4JX6qgcB9ZjYfqAU8sHsvSURERPaGCnU7ufsoYFSZtueB58vZdxqxyynLts8l1hMhIrKNooK8QP02X1SQR2p6jR3vKFJN6A6OIpJ0sQ/mo5IdY4vU9DnJjiASKFobQkRERBJSsSAiIiIJqVgQERGRhDRnQUSkmssvLOqYkZZa6Z8H+YVFkYy01LmVfV7Z+1QsiIhUcxlpqeFjbp5a6eeddWuvCn3GZGdnp/fu3fvwPn36rLr33nu33F75s88+y7zwwgvb33DDDYsGDRr0c0Wfd+zYsQ1nz55d58EHH1y0vX2uvvrq1gMHDszp0qXLptLtw4cPb9OlS5fcnXm+6kDFgoiIJF2dOnUiH3/8cd1IJEI4HPtomjBhQoOsrKzInni++++/X2s+7ATNWRARkaTLzMyMHnzwwZtmzJhRp6Tto48+qnvkkUeuj0ajod/97ncHlLTfddddze69996mubm5KX/4wx/anHHGGe169+7d/sUXX/zF4lD/+9//ap111lmH9urVq/0FF1zQ9ttvv80A6Nevn73//vt1otEoN9xwQ8uTTjrpsH79+tlPP/2UsXdecdWiYkFERAKhd+/eayZNmlQf4OOPP6550EEHbU5LSyvOy8sLffrpp3Vzc3NTotEoU6dObdi/f/+fx4wZ06x9+/abJk2aNP/FF1/85sknn2y2cOHCLYtR5Ofnh6677roDR40a9cPUqVPnXXDBBTl//OMfDyz9nP/5z3/qL1iwoObUqVO/fvjhhxcuXbpUxUI5VCyIiEggnHHGGWs/+uijrKKiIt54440Gffr0WQ1Qs2bNaNeuXde99tpr9T/44IPazZs3z2/RokXh7Nmz67766quNe/fu3b5///6H5uXlpcybNy+z5HwLFizIqF27dtHRRx+9CaBfv35rli5dmrF27drUkn0++uijOieffPLa9PT04v322y/StWvXdXv/lQef5iyIiEgg1K1bN3rQQQdt+uCDD2p/8skndW666aYlEydObABw/vnnr3rkkUeatWjRIv/Xv/71KoBoNBoaPXp0dufOnTcBLF++PNygQYOi8ePHNyjZXvY5iouLKSoq2vI4FAoVR6PRLY9TU1PLHiKoZ0FERAKkd+/ea+67776WZrYpLS1tS/sJJ5ywYdWqVemffvppnb59+64F6Ny58/px48Y1Bli6dGla3759O/zwww9bhiHMLG/9+vWpH3/8cU2Al19+uf5+++1X0LBhwy3VwnHHHZf7zjvv1M/LywutXr069eOPP667t15rVaKeBRGRai6/sChS0cscd/a8GWk795v6GWecsfb2229vPXz48J/KbjvxxBPXrFu3LlyjRo1igBEjRiwdOXJk69NOO61DNBpl+PDhSw4++OD8mTNn1gaoUaNG8d133/39X/7yl/3z8vJS69atG7n//vu/L33Ovn37rp07d27N3r17d2jQoEFh69at83bjJe+zQsXFxcnOkJCZtQGyp02bRsuWLZMdR0T2mOAsJAVVfyGpJUuWcMoppwAc4O6LSm+bO3fuoo4dO65KSrBdEI1GKSgoCA0cOLDt9ddf/2PJsINUvrlz5zbq2LFjm7LtGoYQEZFAW758eVq3bt06dujQYaMKheTQMISIiARa8+bNCz/99NPPk52jOqtwsWBmdwON3X2wmXUCngCygOnA79w9Ymb7A+OA/QAHBrr7BjOrBzwHHAjkABe4+/JKfSUiIiKyR1RoGMLMTgEGl2oaBwx397ZACBgab38YeNjdDyU26Dcq3n47MMPd2xErMv6++9FFRERkb9hhsWBmDYA7gL/GH7cGMt19VnyXZ4DzzSwNOAF4pXR7/Os+xHoWAF4ATo/vLyIiIgFXkZ6Fx4AbgDXxx82BZaW2LwNaAo2A9e4eKdO+zTHx7euBxruVXERERPaKhMWCmV0G/Oju00o1/+KOWEA0QXuiY0REJMmK8/I6Ap0r+0/8vBX25Zdf1jCzzq+++mq9krbFixenDxo06OBevXq1P+200zoMHTr0wBUrVoQhthS1mXUuu4DUQw89tJ+Zdc7Ozk5PdO49pXv37odnZ2enjx07tuHw4cPb7Onn2xt2NMGxP9DMzD4HGgC1gWKgaal9mgFLiU1crGtmqe5eVKod4Kf4MUvMLAzUBbRWuIhIAIRq1Aj/1KJVpZ+3xU8/7tQVdy+++GKj7t27r3nppZcan3vuuWsBrr/++tZ9+/b9uX///qsBxowZ03TkyJGtn3rqqYUAjRo1Knzrrbfql2wHeO+99+rXqlWraEfnlopL2LPg7qe5+2Hu3gm4CZjg7kOAPDM7Lr7bIGCyuxcCM4gVGFva419Pij8mvn1GfH8REREKCwt56623Go4YMeKn7777ruZ3332XAbB69eq0zZs3b/msGjp06MqLL754Zcnjjh075rp7zQ0bNqQALFq0KD0zM7OodLGwvXOX1r1798O3l+3nn39OHTJkyEGnnHJKh969e7d/99136wBMmTKl7llnndXu9NNPb3/ppZcetGrVqu3ernLUqFEte/Xq1f6MM85od9dddzXb2e9Psu3qTZkGAveZ2XygFvBAvP0KYJiZzQO6AzfG20cBx5jZ1/F9rtz1yCIisq+ZPHlyvSZNmhSYWf7xxx+/duzYsY0BrrrqqiUPPfRQy27duh3xhz/8oc3kyZOzTjzxxNyS41JTU4u7dOmyfsqUKVkAr732Wv1evXqtqci5K+quu+5q0apVq/xp06Z9/be//S3773//e4uVK1eG77///pbPPvvsgsmTJ8/r1q3buttvv73c2wwvWrQofdasWVlTp06d9+9///ubH374ocbmzZvLG54PrAp3Ebn7M8SucMDd5wJHl7PPYqBHOe2rgb67mFFERPZxr776asOePXv+DNCnT581N9xwwwE33HDDT717915/4oknzv3vf/9b58MPP6z7wAMPtJw8eXKDZ555ZmHJsX369Fn90ksvNe7Xr9+a999/v/4zzzzz7eOPP958R+dOTU0tPvPMM9tDrAejd+/e7QEef/zxb/fff/8tvd+ff/55nXvuued7gCOOOGLz66+//s2bb76ZlZOTk37RRRcZxG5JXadOnW2GPkq0aNGiID09PXr22Wcf2r1797UjRoz4KTMzM9hrLZShOziKiEhSrVixIjx79uysBQsW1HrppZeaFBcXs2HDhtRXX321/pdfflnrr3/964+9e/de37t37/XXXHPNsu7dux+xcuXKLZ9fPXr0yL311lvbfPnllzWysrIK69WrV7Sjc7/22mv1+/fvv3rKlCnzIDYMUfJ1WeFweJsP9vnz59eIRqOhww47bMOzzz77HcDmzZtDubm55Q5DpKWl8frrr8+fPn16nffeey9rwIABhz777LNuZvmV8f3bG7Q2hIiIJNVLL73UsFOnTrkzZ878YsaMGV9+8MEHX/7mN79Z/vLLLzf+8MMPs8aNG9ewZN+FCxdm1KtXL9KgQYOSy/QJh8N06dJl/ahRo9qcfvrpayp67orm69SpU+7rr7/eAGKFwrBhww456qijNs6bN6+Wu2cA3HPPPc3/8pe/lDsM8emnn2aef/751r1799y//OUvS/bff/+8BQsW1NjZ71MyqWdBRKSaK87Li+zslQsVPW+oxo4/E994442GZZekvvTSS1c+//zzTf71r399c88997R49NFHm2dkZEQbNmxY+PDDD38bDm8bt0+fPquvvvrqg/v06bO2oueeP39+jXbt2uUBzJgx48vt5bvuuuuWjhgxonWvXr3ap6amFt9xxx3ZzZs3L7zpppsWXXXVVQdFo1EaN25c+Pe///378o4/8sgjNx9++OEbe/Xq1aFGjRrRQw45ZFPv3r3X7fAbEyBaolpEAkJLVFemfWmJatl7tES1iIiI7BIVCyIiIpKQigURERFJSMWCiEj1E41Go1XqpkCy58V/Jspdt0nFgohI9fNVTk5OlgoGKRGNRkM5OTlZwFflbdelkyIi1UwkErls+fLlTy5fvvww9EujxESBryKRyGXlbVSxICJSzXTu3HklugW/7ARVlCIiIpKQigURERFJSMWCiIiIJKRiQURERBJSsSAiIiIJqVgQERGRhCp06aSZ3Qb0A4qBf7r7vWZ2KnAvkAm86O43xvftBDwBZAHTgd+5e8TM9gfGAfsBDgx09w2V/HpERESkku2wZ8HMTgROBo4gtobscDPrCDwF/BpoB3Qxs9Pjh4wDhrt7WyAEDI23Pww87O6HElv/dVRlvhARERHZM3ZYLLj7f4GT3D1CrFcgDNQDvnX37Hj7OOB8M2sNZLr7rPjhz8Tb04ATgFdKt1fi6xAREZE9pEJzFty90MxuBeYB04DmwLJSuywDWiZobwSsjxcWpdtFREQk4Co8wdHdbwYaA62AQ8rZJUps2GFn2kVERCTgKjJn4dD4pEXcfRPwKnAS0LTUbs2ApcBP22nPAeqaWWqZdhEREQm4ivQsHAg8YWYZZpZObFLjY4CZ2cHxAuAiYLK7LwbyzOy4+LGD4u2FwAygf+n2ynwhIiIismdUZILjJGAS8BnwCTDT3ccDg4F/E5vH8A1bJy8OBO4zs/lALeCBePsVwDAzmwd0B26svJchIiIie0qF7rMQn69wc5m2aUDHcvadCxxdTvtioMcupRQREZGk0R0cRUREJCEVCyIiIpKQigURERFJSMWCiIiIJKRiQURERBJSsSAiIiIJqVgQERGRhFQsiIiISEIqFkRERCQhFQsiIiKSkIoFERERSUjFgoiIiCSkYkFEREQSUrEgIiIiCalYEBERkYRULIiIiEhCKhZEREQkoXBFdjKzm4EL4g8nuvt1ZnYqcC+QCbzo7jfG9+0EPAFkAdOB37l7xMz2B8YB+wEODHT3DZX5YkRERKTy7bBnIV4U9AR+BXQCOpvZAOAp4NdAO6CLmZ0eP2QcMNzd2wIhYGi8/WHgYXc/FJgDjKrE1yEiIiJ7SEWGIZYB17h7gbsXAvOBtsC37p7t7hFiBcL5ZtYayHT3WfFjn4m3pwEnAK+Ubq+8lyEiIiJ7yg6HIdz965KvzewQoD/wALEiosQyoCXQfDvtjYD18cKidLuIiIgEXIUnOJpZB+Bt4FpgYTm7RIkNO+xMu4iIiARchYoFMzsOmAb82d2fBX4CmpbapRmwNEF7DlDXzFLLtIuIiEjAVWSCYyvgNeAidx8fb/4otskOjhcAFwGT3X0xkBcvLgAGxdsLgRnEhjC2tFfeyxAREZE9pSKXTl4L1ADuNbOStkeBwcC/49smsXXy4kDgCTOrA3xGbH4DwBXAs2Z2I/ADMKAS8ouIiMgeVpEJjlcBV21nc8dy9p8LHF1O+2Kgx07mExERkSSr0E2ZRET2pOLifEKhOcmOsUUsT0ayY4gEhooFEUm6UCiDx349Ltkxtrj89YuTHUEkULQ2hIiIiCSkYkFEREQSUrEgIiIiCalYEBERkYRULIiIiEhCKhZEREQkIRULIiIikpCKBREREUlIxYKIiIgkpGJBREREElKxICIiIgmpWBAREZGEVCyIiIhIQioWREREJCEtUS0iSRcpiARqWehIQYRwut4eRUpU+H+DmdUFZgJnuvsiMzsVuBfIBF509xvj+3UCngCygOnA79w9Ymb7A+OA/QAHBrr7hsp8MSJSNYXTw/R9rU+yY2wx4eyJyY4gEigVKhbMrCuxAqBt/HEm8BRwIvAjMNHMTnf3ycQKgsvcfZaZ/RMYCjwCPAw87O7jzWwUMAr4U2W/IBGpevKLCgL1AZ1fVEBGanqyY4gERkV7FoYCVwL/ij8+GvjW3bMBzGwccL6ZzQMy3X1WfL9ngFvN7EngBODsUu3/RcWCiAAZqekcc/PUZMfYYtatvZIdQSRQKlQsuPtlAGZW0tQcWFZql2VAywTtjYD17h4p0y4iIiIBt6tXQ4TKaYvuQnu1EolEdrzTXhS0PCIiEky7Ot33J6BpqcfNgKUJ2nOAumaW6u5FpdqrlXA4zOOPP57sGFsMGzYs2RFERKQK2NWehY8AM7ODzSwVuAiY7O6LgTwzOy6+36B4eyEwA+hfun03couIiMheskvFgrvnAYOBfwPzgG+AV+KbBwL3mdl8oBbwQLz9CmBYfBJkd+DGXY8tIiIie8tODUO4e5tSX08DOpazz1xiV0uUbV8M9NjphPuQSCQSqK7/SCRCOKwbz4iISGL6pNiLwuEwebe9n+wYW9S4qUeyI4iISBWgtSGkygja1RtByyMisqeoZ0GqDF1NIiKSHOpZEBERkYTUs7AXFRcWBWqeQHFhEaG01GTHEBGRgFOxsBdFiwtJJTgfzkHLsyO6mkREJDn0TrcXpabXAI5KdowtUtPnJDvCTtHVJCIiyaFiQaoMDeOIiCSHioW9qLg4n1AoOL/Nx/JkJDtGhQVt2CRoeURE9pR9qlgI2hhy2TyhUAaP/XpcEhNt6/LXL052hJ2iYRwRkeTQpZMi1URxNDirwgcpi4jsWHB+Da8EQZ8AFymIBOq3+UhBhHD6PvUjIAmEUlIC8/8jSHNPRGTH9EmxF4XTw/R9rU+yY2wx4eyJyY6wTwn6MFiQ5RcWMevWXsmOsUV+YREZmrwqskXVeCcRqQKCfjvq4khwriYpjhQRCm/9MM5IS+WnFq2SmGhbLX76MdkRRAJlnyoWgvRmCL98Q5TdE/SrSYJ+06jYz2IwJoiGwsH5dxSRHdunioUgvRlC1XtDDFq3dVW7miToc2Zk1wX9/4bInrZP/bQF/TfP/KKCQM0TyC8qICM1PdkxKiySH7AJovkRwhn71H8hEZFy7dV3OjO7CLgRSAfuc/d/VOb5iwpTCQfos69snozUdI65eWryApVRdkJZKqEkJSlf2TzhjHCgv38aBtt3qddIqru9ViyYWQvgDqAzkA/MNLP33H1eZT1H0K82CPyM79QIBOmOhEHLsyNBy1smT5B63qra3UNVCEp1tzd7Fk4F3nX31QBm9grQD7htB8elAixfvnyHT1AYLeDhox7bzZiVJ/uH70lL2barY8Ux3ZKU5peazJq5zePCaMEv8iZT2TwFkSivDO2QxETb+n7RD6SHt97XrKiwiNQAXW5XXp5z7puepDTb+s8fT/hF24pQcHq2ipcs2fZxcQGhUHD+b1QkT6n3zOD8UEqVtTeLhebAslKPlwFHV+C4ZgADBw7cE5n2vlqZyU6w1SmnJDuBVFOnvFlOo/5v7CnNgIXJDiFV294sFsr7taEi93ydDXQnVlwUVWoiEZF9VyqxQmF2soNI1bc3i4WfiH3ol2gGLN3RQe6eD3ywp0KJiOzD1KMglWJvFgvvALeYWWNgI3AeEJw72IiIiEi59tqqk+7+E3AD8B7wOfC8u3+8t55fREREdk2ouLg42RlEREQkwPZaz4KIiIhUTSoWREREJCEVCyIiIpKQigURERFJSMWCiIiIJKRiQURERBJSsSCVxsxCZtYg2TmkejGz25OdYUfMLD3+98Fm1sfM9N4rVYrusxBwZtYReBZoBfwH+KO758a3feruRyYxWyvgLmA18ATwBpAJ5AD93H1+srJVBWaWBdxK/N/W3ceV2va4uyf1DqdmFgYGA2uAt4BHgcOJ3X79TyU/h8lmZnOBTu4eyDczM7sJOBi4EZgFzAOy3X1oUoOJ7IS9ebvnQDKzp4Htvsm4+6V7MU55Hgb+CHwB/AV4z8x6uPsGyl+ca296BngJaA28Dwxw96lmdhKx3CclLxqY2aBE29197N7Ksh1PA18CM4A/m9kJpQqEo5IXa4vHiRV/+wGjgInAHcAFxAqHoCwF+zPwjZl9CmwuaQzA/90SfYHjiP0/Hufu15nZnCRnEtkp1b5YIPZGfT8wAshLbpRy1XT39+JfX2Fm9wATzKxXMkPFNXT3x+Jdqpe4+1QAd3/PzMYkORvAyUA/YgVN2cKqGEh2sXCAu58LYGaTgIlmNsbdryH5hSBAF3c/3MxqAYvd/YZ4+y1m9lkyg5XxbLID7ECqu+eb2ZnAjfH/L7WSHUpkZ1T7YsHdnzKztsTeuP+c7DzlyDWz04Ep7l7s7tea2XPAv4GaSc620cxOc/e3zax9SaOZnU1ssbCkcvfB8TkUH7j7U8nOUx4za+ruy919s5mdA0w3s+tJ0Nu1F0XNrJG7rzKzi0sazawlAZrv5O7PmlkboAMwFWjl7tnJTbWNaWb2FbAJmA78F5iQ3EgiOycw/+GT7GZgWrJDbMflwPXAb0q1DQK+Bw5MSqKthgIjzSzF3dcBmNn5wMj4tiC4HAjqpMtbgE/MrC9A/HvYm9iKrEckMVeJW4DPzCzV3acAmNlpwCfEhiUCwcz6E5sv83di/9b/K13cBMAdwBnAse4eBYYTG6YTqTI0wbEcJb/tJTvH9pTkK/mtL9l5Sgv69y5ozKwOkObuq0u1pQB93f21pAXbmqWmu28q9bgekFI6b7LF5yqcCEx391+ZWTPgHXfvkORcrYgNJ00CTmfr0FIYmOTuhyYrm8jOUs9C+SYlO8AOTAIIWqEQF+jvnZm9mewMpbl7bplC4U13jwahUAAoXSjEjQtSoRBXVPrKDHdfBkSTmKfErcSGHA5h6/DDf4kNlUxOYi6RnVbt5yxsRxAmlyUS5HxBzgbQPNkBdkD5dt7XZvYHIM3MOgFXAJ8nNRFbr8Ywsz+5++j416GgXuIpkoh6Fsq3NNkBdiDI+YKcDYJfzCjfzrsSaEHsssmngPXECoag+NjMPox/3dbMvjezbklNJLKTNGdB9nlmdpS7z4l/3SzeTR0Yyrdvi8+pGOTuX8UfHwr8y927JDeZSMWpWIgzs6HEZi03jDeFgGJ3T01eqq2CnC/I2QDMbDLQFngPeBN4q5yx+KRRvt1TBX7+5rl7+zJtn7t7pyRFEtlpKhbizCwbONPdv052lvIEOV+Qs5UwsxrE7ijZGzgLcHc/PbmptlK+XRf0nz8zexX4FvhXvOlCoK27X5C8VCI7RxMct1oZ1DebuCDnC3I2zKwxsUvregDdia1lEZi8yrfbAv3zB/wWuB14ASgkdmVEUO5DIlIh1b5nodT6AWcCGcDrQKRke7LXDwhyviBnK83MosAKYrf1ftzd1yQ30baUb9dUlZ+/sswsROyOsd8nO4tIRalnYetiRxvjf7qX2haE9QOCnC/I2Uoz4BRied8zs3nAe+7+RHJjbaF8u6ZK/PyZ2XBicypKrwexCDgoKYFEdkG171mQ6sPMjgBOA34H4O6HJDfRtpRv15SsT1Km7Vx3fzVZmUqLz6k4mVjBcD2x4ZzT3D0oq3aK7JB6FuLM7Fug9OzpYmLXbc8HrnX3xUkJFhfkfEHOBmBm44ktEfwNsTtMnununsxMpSnfromvCZEB3GZmN5XaFCb2oRyIYoHYnIpsM/sCONzdn4nfREqkylCxsNVkYoszlaxOOBDoQmyBmn8CpyYpV4kg5wtyNogtUX0ZsUvqUt19bXLj/ILy7Zq6QDegNrHf1ksUATeUd0CSbDSzk4AvgLPNbDZQP8mZRHaK7uC41fHufr+7r4//eQQ4wt3/QzBWLQxyviBnA5gLvEtsnPh7M/ssvix5UCjfLojPmXgGqEdspc7zgFbEbnj0YvKS/cJwoC8whdi9IL4BHkxqIpGdpGJhqyIz61XyIP51gZk1AdKSF2uLIOcLcjaAR4G/uXtDd28A3Ak8nuRMpSnfLjCzk4HniX0YdyPWu/Aa8IKZ9UhasDgzezf+5Tnu/sf4AmHnuXs9d78/mdlEdpaGIbYaAjxjZs8R6279FhgMDAPuSWKuEkHOF+RsAI3c/ZWSB+7+kpndmMxAZSjfrrkZ6OPun5dq+8zMZgH3ASckJdVWbczsduDS+LLj23D325KQSWSXqFiIi9+3/Sgzq09sydv18U1/SWKsLYKcL8jZ4vLN7Eh3/xTAzDoDgbldMcq3q+qWKRQAcPdPzCwIw1/nEbsHRIhgLsAlUmHVvlgws8fdfZiZvUdsFn9JOwDufnKyssVzBDZfkLOVcTXwbzNbTexNuwGxW+4GxdUo366obWZhd4+UbjSzMAF4b3P3z4j1dMxx98nJziOyO5L+HyoAHov/fUsyQyQQ5HxBzraFu8+KT8hrS2yejrt7QZJjbaF8u2wqMBq4pqTBzFKJDUFMTFaoUlked/dhwHVmNqLs9gAV0yI7pJsylWJm7YBGlOoydPfpyUu0rSDnC2I2M3uaUj0eZbn7pXsxzi8o3+4xs1rELs/dH5hD7Jefo4itW3Guu+cnMR5m1jk+JHJiedvd/b97O5PIrlLPQpyZPQacASxk6xtkMbE7ryVdkPMFONv7SX7+HXk/2QF24P1kB0jE3TcCJ8c/jLsQ+5m7390/SG6yGHf/JP5lP3cfXnqbmT0LqFiQKkM9C3FmthBoF5Du1V8Icr4gZyvLzIa5e9Iv+9se5dt3mNmTwIHEejvmlNoUBuq5+xFJCSayC9SzsNUPQCYQ1A+8IOcLcrayfkcA7hGQgPLtO24H2gB/B24t1R4hdit0kSqj2hcLpcZlw8BcM5vOtsvcJntcNrD5gpwtgaBfwqZ8+wh3XwQsMrMVmp8gVV21LxbYOi4b1P/M78f/DmK+9+N/BzHb9jya7AA7oHz7nhpm1srdf0x2EJFdpTkL5Qj6uGyQ8wUpm5k9ANzs7muSnaU8ylc9mNk3wCHASmKrsQLg7gcmLZTITtLaEOX7XbID7ECQ8wUp2yBglpmdm+wg26F81cPZwAjgEWJzF24HPkxmIJGdpWKhfEEflw1yviBlywbOAa4ys4/MrL+ZZSY7VCnKVz2MAc4ChgK9gNuAGklNJLKTVCyUL+jjskHOF6Rsxe4+z91PBG4gdq/+bDObbmbPJzkbKF91YcTuOfIf4G/A0UCLpCYS2UnVvlgwswfiCyBt4e6PbW//vS3I+YKcLa703STfcfcLgFbAH4m9cSeb8lUPK929GPgGOMLdlwIZSc4kslN0NURsXLaXmY1091eTHaYcQc4X5GwAD5VtcPdC4JP4n2RTvurhKzN7kNichefMrDmQluRMIjtFxUJsXHYg8IiZ/Qm4F5jg7psTH7bXBDlfkLMBfFHyhZmdQuyW1IXAf9z9o6Sl2kr5qoffA93cfZ6Z3QycAlyU5EwiO6XaD0MQ/HHZIOcLcjaIr4ppZlcC9wM/AiuAx8zsD0nMVUL5qgF3L3L3GfGvJ7j7Ve7+VbJziewM9SyUGZcF3jGzNOAIYvd1T7Yg5wtyttKGAj3c/WfYcs/+2ZTTzZ4kyicigaZiIfjjskHOF+RsAGlmlkLsZjgbS7UXANHkRNqG8olIlaBiIfjjskHOF+RsADnEus6LiV3SOdjMTiZ2+drLyQwWp3wiUiVozkLwx2WDnC/I2XD3k929BbEJZSX3f8gndgvjm5OXLEb5RKSqUM/CVkEflw1yviBng9hlai3NrA8wz92Ddqtd5RORQFOxEPxx2SDnC3I2zGw/4BXgMOBbYt3pZmb/Ay5y97VJjKd8IlJlaBhi67hse+JdrfFx2Q8JxrhskPMFORvAg8AHQBN37+ruxwBNgLnEhk2STflEpEqo9j0L7n4yxH5dAkpuXVwyLjsxacHigpwvyNnijnD3/qUb3L3AzK4HPk9OpG0on4hUCepZ2Kr0uOzSgHzYlRbkfEHNlldeY/w+/UkfJkH5RKSKqPY9C0Eflw1yviBniyvexW17i/KJSJVQ7YsFto7LnhK/oRBmlg7cSmxcdnDSksUEOV+QswF0MLPvy2kPAc32dphyKJ+IVAkqFoI/LhvkfEHOBtAWaAikErtiA+Ak4OtSj5NJ+USkStCcheCPywY5X5CzATQAJgIN3X2xuy8GDgFeA7KSGSxO+USkSlCxEPxx2SDnC3I2gHuAAe4+paTB3W8ALiW2nHayKZ+IVAkahgj+uGyQ8wU5G0B9d3+/bKO7TzWz0UnIU5byiUiVoGIh+OOyQc4X5GwQv8Oku28zJBK/62R6kjKVpnwiUiVoGCL447JBzhfkbAD/Bcpb8OhGYM5ezlIe5RORKiFUXByEoeXkMbNpwF/KdreaWS9ghLufmpRgW3MENl+Qs8Vz1AEmERsSmU1seORIYr0efd19dRLjKZ+IVBkqFsw+dfcjt7Ptc3fvtJcjlc0Q2HxBzlYqR4jY0MiviF2hMcfdZyQ31VbKJyJVgeYsBH9cNsj5gpwN2HIZ57vxP4GjfCJSFWjOQvDHZYOcL8jZRESkkmgYIuDjskHOF+RsIiJSeap9sQDBH5cNcr4gZxMRkcqhYkFEREQS0pwFERERSUjFgoiIiCSkYkFEREQS0n0WpMows2LgK6CI2KqWNYH1wO/dfU58n1rArcBZQEF8vzeA2919c3yf94GH3P2VBM91FPBnd++3E/kaATnuHjKzvsCp7v5/CfZ/Ehjv7u+U89yvuHubij53mePbAF+5e+1dOV5EpCwVC1LVnOTuq0oemNm1wIPAsWYWBt4B/gf8yt03mVlN4E5gqpmd7O6RijxJvPiocKFQzvETgAk72OeyXT2/iMjepGJBqqx4cbA/UHI/h/OBFHf/fyX7xAuGq4HPgHOAlyt47h7Eeh8OM7NniPVgHA60Ar4BLnT3DWZ2LnAHsInYvSZKjh9MrNj4f8BMoLm7F5hZKrAY6Ak8HH+OV8zs98AfgXXAl6XOcwvQyN3/UPaxmR0D/A3IIHavi7fd/bcVeX0iIjtDcxakqnnPzOaa2VJgQbxtSPzvbsD0sgfEb1k8DTh+N563M9AbaAc0B843sybAU8B57t6ZWBFQ9rkXEFuyu2+8qSewyN3nlexjZp2AW4AT3L0LseGTirgKuMnduwLtgb5m1nnnX5qISGIqFqSqOcndOwJ9iM1ZmOnuK0ttT9vOcRnE5i/sqinunu/uhcR+829ArPj4stQH/2PbOfYJYHD86yHAk2W2nwK85e7L448fr2CmS4B6ZnY9sV6KmoDmKYhIpVOxIFWSu39GrNv+yfiEPoAPgRPiC1ltEX98ArHhgF21udTXxcRubV3yd4ntzYd4BehqZu2AE4GXymxPdJ6y20ov0DUDOIPYsMhtwJIy+4qIVAoVC1JlufsLxCYz3h9vegXYCNxvZpkA8b8fBDYA/6nkCDOADmbWMf548HZy5gHjgWeAf7v7pjK7vA30NLOW5ZwnB+hsZqH4lR49AcysPnAU8Cd3fxVoARwMpO7maxIR+QUVC1LV/QE43cx6xa906EmsMPjEzL4CPo0/Pi0+hFBp3D0HuAh4zsw+BQ5IsPsTwNH8cggCd/8SuA6YZmZzgBqlNj9HrGD4ltiiXf+LH7OG2FUen8aPGUmsZ+Xg3XxZIiK/oLUhREREJCFdOinVlpm9CNh2Nvd3d9+beUREgko9CyIiIpKQ5iyIiIhIQioWREREJCEVCyIiIpKQigURERFJ6P8Dl38UQaB1u/8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# https://scimap.xyz/All%20Functions/C.%20Plotting/sm.pl.stacked_barplot/\n",
    "sm.pl.stacked_barplot (adata,\n",
    "                       x_axis='ROI_individual',\n",
    "                       y_axis='phenotype',\n",
    "                       method='absolute')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0404b8ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the number of cells normalized to 100% \n",
    "sm.pl.stacked_barplot (adata,\n",
    "                       x_axis='ROI_individual',\n",
    "                       y_axis='phenotype',\n",
    "                       method='percent')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00776b3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# specify the elements to be in the plot\n",
    "x_axis_elements = ['CD57-low-1', 'CD57-low-2', 'CD57-low-3', 'CD57-high-2', 'CD57-high-1', 'CD57-high-3']\n",
    "y_axis_elements = ['ASMA+ cells', 'Myeloid', 'NK cells', 'Neutrophils', 'Other Immune cells', 'Treg', 'Tumor']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28dfe1a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# replot\n",
    "sm.pl.stacked_barplot (adata,\n",
    "                       x_axis='ROI_individual',\n",
    "                       y_axis='phenotype',\n",
    "                       method='percent',\n",
    "                       subset_xaxis=x_axis_elements,\n",
    "                       subset_yaxis=y_axis_elements)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68337884",
   "metadata": {},
   "outputs": [],
   "source": [
    "# quiet a number of parameters to play around:\n",
    "sm.pl.stacked_barplot (adata, \n",
    "                x_axis='ROI_individual', y_axis='phenotype', \n",
    "                subset_xaxis=x_axis_elements, subset_yaxis=y_axis_elements, \n",
    "                order_xaxis=None, order_yaxis=None, \n",
    "                method='percent', plot_tool='plotly', \n",
    "                matplotlib_cmap=None, \n",
    "                matplotlib_bbox_to_anchor=(1, 1.02), \n",
    "                matplotlib_legend_loc=2, \n",
    "                return_data=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3430ccd2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b4a9013f",
   "metadata": {},
   "source": [
    "### Calculate the fold change in cell types between the different ROI's"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebca1965",
   "metadata": {},
   "outputs": [],
   "source": [
    "adata = sm.tl.foldchange (adata, \n",
    "                          from_group=['CD57-low-1', 'CD57-low-2', 'CD57-low-3'], \n",
    "                          to_group=None, \n",
    "                          imageid='ROI_individual', \n",
    "                          phenotype='phenotype',\n",
    "                          normalize=True, \n",
    "                          subset_phenotype=None, \n",
    "                          label='foldchange')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecef1c39",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Heatmap of foldchnage  \n",
    "sm.pl.foldchange (adata, label='foldchange', method='heatmap',\n",
    "                     p_val=0.05, nonsig_color='grey',\n",
    "                     cmap = 'vlag', log=True, center=0, linecolor='black',linewidths=0.7,\n",
    "                     vmin=-5, vmax=5, row_cluster=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c0d912f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parallel_coordinates plot of the foldchanges\n",
    "sm.pl.foldchange (adata, label='foldchange', \n",
    "                  subset_xaxis = ['ASMA+ cells', 'NK cells', 'Neutrophils', 'Treg', 'Tumor'],\n",
    "                log=True, method='parallel_coordinates', invert_axis=True,\n",
    "                parallel_coordinates_color=['black','blue','green','red','#000000'],\n",
    "                matplotlib_bbox_to_anchor=(1.04,1),\n",
    "                matplotlib_legend_loc='upper left',\n",
    "                xticks_rotation=90,\n",
    "                return_data = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d29c09d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save adata\n",
    "adata.write(str(common_path) + 'may2022_tutorial.h5ad')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
